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ABSTRACT 

Motivated by the recent results of lLesur fc Qgilviel ()2010T ) on the transport prop- 
erties of incompressible convection in protoplanetary discs, in this paper we study the 
role of compressibility and hence of another basic mode - spiral density waves - in 
convective instability in discs. We analyse the linear dynamics of non-axisymmetric 
convection and spiral density waves in a Keplerian disc with superadiabatic verti- 
cal stratification using the local shearing box approach. It is demonstrated that the 
shear associated with Keplerian differential rotation introduces a novel phenomenon, 
it causes these two perturbation modes to become coupled: during evolution the con- 
vective mode generates (trailing) spiral density waves and can therefore be regarded 
as a new source of spiral density waves in discs. The wave generation process studied 
here owes its existence solely to shear of the disc's differential rotation, and is a special 
manifestation of a more general linear mode coupling phenomena universally taking 
place in flows with an inhomogeneous velocity profile. We quantify the efficiency of 
spiral density wave generation by convection as a function of azimuthal and vertical 
wavenumbers of these modes and find that it is maximal and most powerful when 
both these length-scales are comparable to the disc scale height. We also show that 
unlike the convective mode, which tends to transport angular momentum inwards in 
the linear regime, the spiral density waves transport angular momentum outwards. 
Based on these findings, we suggest that in the non-linear regime spiral density waves 
generated by convection may play a role in enhancing the transport of angular mo- 
mentum due the convective mode alone, which is actually being changed to outward 
by non-linearity, as indicated by above-mentioned recent developments. 
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- (stars:)planetary systems: protoplanetary discs - turbulence 



1 INTRODUCTION 



The main agent responsible for the anomalous outward 
transport of angular momentum in protoplanetary accre- 
tion discs is widely recognised to be turbulence due to 
magnetorota tional instability (MRI. lBalbus fc Hawlevll 19981 : 
Balbus 2003). However, in order for the MRI to operate, 



magnetic field are effectively coupled 


Blaes & Balbus! 1994 


Sano et all l200d; 


Sano & Stone 2002 


; Salmeron & Wardld 


20031: Desch 2004 


). In cold and dense interiors of proto- 



planetary discs, the necessary level of ionisation is not typ- 
ically reached and, as a result, a large magnetically inac- 
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tive region -'dead zone' - develops in the disc (se e e.g., 
lGammie![l99d : IStone et al1l200d : iFromang et al.ll20o3 ). Sev- 
eral mechanisms have been put forward to explain (non- 
magnetic) transport in the dead zone: penetration/diffusion 
of turbulent fluctuatio ns into the dead zone due to the MRI 
in the surface layers (Fleming fc Stond 120031: | Turner et al.l 



120071 : lOishi et all 120071: lOishi fc Mac Lowil2009l ). transport 
due to self-gravity l|Armitage et al.l l200lF ) and due to vari- 
ous hydrodynamic candidates, such as vortices, spiral den- 
sity waves, Rossby wave and baro c linic instabilities (e.g., 
Lovelace etafl Il999l: [Li et alj l200ll ; iKlahr fc Bodenheimerl 
20031: iJohnson fc Gammid l2005bl: iPetersen et all 120071 : 



iJohnson fc Gammidl2006l : iLesur fc Paparoizoull2010f ). 

Recent developments indicate that yet another mecha- 
nism - thermal convection in the vertical direction - can also 
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provide outward transport of angular momentum. Vertical 
convection was first suggested as a source of angular momen- 
tum transport in protoplanetary discs long before the MRI 
llCameronll 19781; iLin fc Papaloizoulll980l . Il985l ; iRuden fc Linl 
ll986l : lRuden fc Porlacklll99lf ). In these earlier studies, it was 
pointed out that the temperature dependence of the opacity 
in protoplanetary discs can lead to convect ive instability for 
a wide range of disc temperatures (see also|Rankov 20071). A 
corresponding effective viscosity was estimated based on a 
phenomenological mixing-length prescription and using this, 
evolutionary models of discs - with outward angular mo- 
mentum transport driven solely by co nvective turbulence - 
were constructed. IRuden et al.l (|l988l ) carried out a linear 
analysis of axisymmetric convective modes in discs and, al- 
though linear axisymmetric modes do not produce torques 
themselves to transport angular momentum, these authors 
also estimated - from the radial wavelengths and growth 
rates of the most unstable axisymmetric modes - a Shakura- 
Sunyaev a ~ 10 -3 — 10 -2 , which might be in the case of 
non-linear axisymmetric or non-axisymmetric vertical con- 
vection. However, the non-linear development of these ax- 
isymmetric (two-dimensional) convective modes turned out 
to lead to inward (i.e . , towards the central star) transport o f 
angular momentum (|Klev et al 1 ll993l ; iRiidiger et al"1l2002t] , 
which is obviously not what is required for the accretion of 
matter from the disc onto the central star. 

To investigate the problem of convective transport more 
fully, subsequent studies considered the dynamics of non- 
axisymmetric convection in discs with Keplerian differential 
rotation, or shear. In this case, shear-induced effects (see 
below) come into play for non-axisymmetric perturbations 
and hence their dynamics is rich er than that of ax isym- 
metric ones. The linear studies of iKorvcans <|1992| ) and 
iRvu fc Goodman! (|l992l . hereafter RG92) in the shearing box 
framework showed that growing non-axisymmetric shear- 
ing waves of the convective mode in the trailing phase lead 
predomin antly to inward transport of angular momentum. 
However, iLin et all l|l993l) argued that this can be changed 
if the radial stratification, which is absent in the shearing 
box, is taken into account. A stratified background can sup- 
port radially localised packets of the convective mode, which 
are able to provide an outward flux of angular momentum. 
Following more gener al non-linear three-dime n sional simula - 
tion s of convection byfCabot fc Pollackl (|l992l ); ICabotl (|l996l ) 
and lStone fc Balbusl (| 19961 ). however, did not yield outward 
angular momentum transport either; the time-averaged a- 
parameter turned out to be small and negative, implying 
inward transport. This led to vertical convection being re- 
garded as an inefficient mechanism for driving angular mo- 
mentum transport and disc's secular evolution. However, 
iKlahr et all (| 19991 ). using global disc simulations, demon- 
strated that convective flow can in fact be non-axisymmetric 
(i.e., have azimuthal structure) and pointed out that the 
small azimuthal extent of the computational domains and 
low Reynolds numbers in the simulations of ICabotl (| 19961 ) 
and IStone fc Balbusl §996) had a tendency to wipe out all 
azimuthal variations resulting in the near axisymmetry of 
convection, which as noted above, is characterised by inward 
transport of angular momentum. 

Since high-resolution simulations with improved numer- 
ical techniques are becoming increas ingly affordable, interest 
in convective transport is reviving. iLesur fc Ogilviei (|2010l ) 



performed incompressible simulations of convective insta- 
bility in the shearing box and found, contrary to previous 
results, that at high enough Rayleigh (Reynolds) numbers 
( Si 10 6 ), the sign of the corresponding a is reversed to 
positive, though it still remains small. This was demon- 
strated to be a consequence of the non-axisymmetric struc- 
ture of convection, which is established at large Rayleigh 
numbers. This plays a central role in that it gives rise to 
an appreciable pressure-strain correlation tensor that, in 
turn, changes the direction of transport from inwards to 
outwards. Past simulations mentioned above, being at low 
Rayleigh (Reynolds) numbers, had thus been unable to cap- 
ture the non-a xisymmetry o f conv ective flow. Analogous 
simulations by Ka pyla, et al.l |20ld ) also showed that due 
to non-axisymmetry, convection can transport angular mo- 
mentum outwards. However, the specific parameter range 
considered by these au thors is somewhat different to that of 
ILesur fc Ogilviei (|2010l ). Although the findings of both these 
studies seem encouraging, they should be viewed as tenta- 
tive and requiring further corroboration. 

All the above-mentioned studies consider perturbation 
dynamics in a sheared environment, as the Keplerian dif- 
ferential rotation of protoplanetary discs is characterised by 
a strong radial shear of the azimuthal velocity. It is well 
known from fluid dynamical studies that operators govern- 
ing the linear dynamics of perturbations in flows with non- 
uniform kinematics are non-normal due to the shear of the 
mean velocity profile, resulting in a numbe r of linear tran- 
sient , or finite-time phenomena ( see e.g., iTrefethen et al.l 
1 19931 ; ISchmid fc Henningsonl l200ll . for an introduction to 
the subject). The standard modal approach (i.e., spec- 
tral expansion of perturbations in time and examination 
of eigenfrequencies), com monly employed in hydrodynam- 
ics (|Drazin fc Reidlll98ll ). describes perturbation behaviour 
(stability) only at asymptotically large times; in fact it is un- 
able to capture the transient dynamics of perturbations at 
intermediate times. In other words, predictions of the modal 
analysis concerning flow stability are really relevant only to 
the asymptotic fate of a flow. Accordingly, new mathemati- 
cal methods - broadly known as the non-modal approach 
- had been developed that allow for the non-normality- 
induc ed perturbation dynamics (e.g.. lSchmid fc Henningsonl 
2001). The non-modal approach leads to an initial value 
problem, enabling us to trace a full temporal evolution of 
perturbations and in that respect it is advantageous over 
the modal approach. 

One of the most important transient phenomenon, 
which is a consequence of the non-normality, is the coupling 
and energy exchange among different perturbation modes 
with each other and with the mean flow occurring during a 
finite time interval in the linear theory. The linear mode 
coupling is a gene ral phenomenon intrinsic to inhomoge- 
neous/shear flows (|Chagelishvili et al.lll997l ) and plays an 
important role in the subsequent non-linear development of 
perturbations and largely defines the characteristics of the 
resulting non- linear state (turbulence). Because of strong 
Keplerian shear, similar non- normality /shear-induced tran- 
sient processes should inevitably take place in discs as well. 
So, for a proper understanding of energy exchange processes 
in the non-linear regime and ultimately of disc turbulence 
phenomenon, it is necessary to first analyse in the linear 
regime all possible couplings and energy exchange channels 
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among perturbation modes existing in disc flows. As pointed 
out above, the modal approach, which has been as widely 
employed in studying perturbation dynamics in discs (e.g. , 



Narayan et al. 1987; Lin et al. 199 01; iLubow fc Ogilvielll99a : 



Ogilvie fc Lubowll 19991 ; iLi et al.ll2003r 7 as in fluid mechanics 



does not account for a finite-time aspect of perturbation dy- 
namics originating from the shear of disc flow. For this rea- 
son, when analysing temporal evolution of perturbations in 
discs, a different technique, a special type of the non-modal 
approach - th e method of shearing waves - is of t en employed 
(e.g RG92; iGoldreich fc Lvnden-Belll Il965l; iBodo et al ' 



2005 



§; Johnson fc Ga mmic 20 05al ; iHeinemann fe Papaloizou 



2009a, hereafter HP09a), which we also adopt in this paper. 

A special manifestation of the shear-induced lin- 
ear mode coupling phenomena in discs is the genera- 
tion of spiral density and inertia-gravity waves by vor- 
tices, which has already been extensively studied tak- 
ing into account the presence of other factors spe- 
cific to discs: radial and vertical stratifications, self- 
gravity, MRI turbulence, etc iDayis et al.ll2o"oo1; |Pavis!l2002l; 
Tcvzadzc et al. 2003; Bodo ct al. 2005; Johnson & Gammic 
2005al lbl; IBodo et all 120071 ; IMamatsashvili fc Chagelishvilil 
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20071; ITevzadze et alj 120081; IMamatsashvili fc Rice! l200g[ 



Heinemann fc Papaloizoull2009al lbl: iTevzadze et alJl201Ch . "in 



particular, it was shown that an efficient generation of spi- 
ral density waves (SDWs) by vortices occurs when the char- 
acteristic horizontal length-scale of these two perturbation 
types is of the order of the disc thickness, which, in turn, 
implies that compressibility effects are important at such 
length-scales. The role SDWs play in disc dynamics cannot 
be overestimated. In particular, it was suggested that gravi- 
tational forces due to stochastic density perturbations asso- 
ciated with SDWs in a turbulent disc fl ow, may be important 
in th e migration of low-mass planets (|Nelson fc Papaloizoul 
120041 ; iNelsonl l2005h . More significantly, SDWs are able 
to enhance angular mom e ntum transport rate ( see e.g., 
I Johnson fc Gammie|[2005b1 ; lOishi fc Mac LowU2009h . SDWs 
also play a role in the dynamics and angular momentum 
transport in a dead zone dFleming fc Sto ne 2003; Oishi et al.l 
120071 ; lOishi fc Mac Low! 120091 '). So. it is crucial to identify 
and analyse the generation mechanisms of SDWs. 

Our goal in this paper is to investigate yet another man- 
ifestation of shear-induced linear mode coupling phenomena 
in discs - generation of SDWs by convection. In a disc with 
superadiabatic vertical stratification, together with the con- 
vective mode due to a negative entropy gradient, there also 
exists a SDW mode due to compressibility. As we will show, 
Keplerian shear of the disc's differential rotation causes the 
coupling of these two modes in a way similar to that of cou- 
pling of SDWs to vortices mentioned above. Since normally 
SDWs contribute to the angular momentum transport pro- 
cess, this coupling can in turn have important implications 
for convective transport as well. So, convection, like vortices, 
can be regarded as a new source of SDWs. However, we omit 
the vortical mode in this paper, because in the linear theory 
presented here, an exponentially growing convective mode 
is more powerful and dominates the algebraic growth of the 
vortical mode (unless the disc is self-gravitating). Follow- 
ing RG92, we treat a stratified disc in the local shearing 
box approximation. We examine in detail how the gener- 
ation of the SDW mode by the convective mode occurs 
during evolution and quantify its efficiency as a function 



of azimuthal and vertical wavelengths of these modes. We 
demonstrate that the efficiency of wave excitation is max- 
imal and powerful when both these wavelengths are com- 
parable to the disc scale height, that is, for wavelengths 
at which compressibility is important. Indeed, the above- 
mentioned numerical simulations indicate that convective 
motions (rolls) actually extend over the entire scale height 
of the di sc and are characterise d by moderate Mach numbers 
(see e.g. JStone fc Balbusll 19961 ). implying that compressibil- 
ity and, hence SDWs, also actively participate in dynami- 
cal processes in convectively unstable discs and their effects 
must be understood properly. Finally, we would like to note 
that SDW generation can also be seen in the linear analysis 
of RG92 and also in a related linear study of the growth of 
non-axisymmetric shearing wave s of the convective mode by 
iBrandenburg fc Dintransl (I2006T ). though these authors do 
not specifically focus on and characterise this phenomenon. 

The paper is organised as follows. Physical model and 
basic equations are introduced in Section 2. Shear-induced 
linear coupling of SDWs and convection is investigated in 
detail in Section 3. Angular momentum transport by these 
two modes is analysed in Section 4. Summary and discussion 
are presented in Section 5. 



2 PHYSICAL MODEL AND EQUATIONS 

To study the dynamics of perturbation modes existing 
in compressible and stratified gaseous discs with Keple- 
rian rotation, we adopt a loca l shearing box approach of 
IGoldreich fc Lvnden-Belj (|l965h - In the shearing box model, 
disc dynamics is studied in a local Cartesian reference frame 
co-rotating with the disc's angular velocity at some fiducial 
radius ro from the central star, so curvature effects due to 
cylindrical geometry of the disc are ignored. In this coor- 
dinate frame, the unperturbed differential rotation of the 
disc manifests itself as a parallel azimuthal flow with a lin- 
ear velocity shear in the radial direction. The Coriolis force 
is included to take into account the effects of the coordi- 
nate frame rotation. The vertical component of the gravity 
force of the central star is also present, but self-gravity of 
the disc is neglected. As a result, we can write down the 
three-dimensional shearing box equations of motion 

^ + 20z x u + V ( -g«V + ifiV | + -Vp = (1) 
at \ 2 J p 

and continuity 



dp 

~dt 



+ pV ■ u = 0. 



(2) 



Here u = (u x ,u y ,u z ) is the velocity of gas in the local frame; 
p and p are, respectively, the gas density and pressure; fi is 
the angular velocity of the local rotating reference frame, 
which is equal to the disc's angular velocity at ro: Q(ro); 
x,y,z are, respectively, the radial, azimuthal and vertical 
coordinates; z is the unit vector along the vertical direction; 
d/dt = d/dt + (u ■ V) is the total time derivative operator. 
The shear parameter q = 1.5 for the Keplerian differential 
rotation considered in this paper. 

In the absence of heat exchange, the specific entropy s = 
c„ln (p/p 7 ) of fluid elements is conserved along streamlines 



(3) 
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where the adiabatic index 7 = c p /c v is the ratio of specific 
heats c p and c„ at constant pressure and constant volume, 
respectively (for convenience, below we replace s/c v — > s). 
It follows from equations (1-3) that Ertel's theorem of po- 
tential vorticity conservation holds 



d (V x u + 2fiz) • Vs 
dt p 



= 0, 



(4) 



which we will make use of later in classifying perturbation 
modes in the disc. 



is the vertical stratification scale height of the disc, p m = 
po(0) and p m = po(0) are the midplane values of the equilib- 
rium density and pressure. As mentioned above, g is some 
height-averaged value of acceleration and we choose g 2 to be 
equal to the density- weighted average of (Q 2 z) 2 , where as a 
weight function we use an exact expression for the equilib- 
rium density, p m exp(— ^yQ 2 z 2 /2c 2 ), which is obtained from 
equations (5) and (6) with a linearly increasing gravitational 
acceleration Q 2 z (RG92). This relates to the angular veloc- 
ity through 



7ff 

„2 



2.1 The equilibrium disc model 

Equations (1-3) have an equilibrium solution that is station- 
ary and axisymmetric. In this unperturbed state, the veloc- 
ity field of Keplerian rotation represents, as noted above, 
a parallel azimuthal flow, uo, with a constant radial shear 
-qQ: 

u-xo = u z0 = 0, Uyo = —qQx. 

In the shearing box model, equilibrium density, po, and pres- 
sure, po, depend only on the vertical z— coordinate and sat- 
isfy the hydrostatic relation 



9o = 



(5) 



1 dp n 2 
= II z. 

Po dz 

Following RG92 an d lBrandenburg fc Dintransl l|2006tl , we as- 
sume the unperturbed disc to be vertically isothermal with 
a constant adiabatic sound speed 



(6) 



For simplicity, we use th e thin dis c approximation with 
constant gravity (see e.g., IShulll974l ; iTevzadze et al.ll2003l . 
RG92). In other words, the vertically varying gravitational 
acceleration ft 2 z is replaced with its some height-averaged 
constant value g > 0, but because the former has opposite 
signs on either side of the disc midplane z = 0, we write 



go = ft z -t sign(z)g. 



(7) 



Although this simplification is highly idealised, it greatly 
facilitates the mathematical treatment of the problem al- 
lowing us to make Fourier analysis of perturbations in the 
vertical direction too in addition to horizontal decomposi- 
tion into shearing waves (see below). This, in turn, makes it 
possible to clearly understand the key effects of shear on the 
dynamics of perturbation modes with various vertical and 
horizontal length-scales in a stratified disc and avoids the 
need to solve complex partial differential equations with re- 
spect to time t and the vertical coordinate z. In this case of 
constant g, substituting replacement (7) into equation (5) 
and using expression (6) for the vertically constant sound 
speed, we find the distribution of the equilibrium density 
and pressure with z: 



where 



po(z) _ po(z) 



H 



exp 



75 



l£l 
H 



(8) 



The Brunt- Vaisala, frequency squared is defined as 

N ^9o(ldpo_dpo\ = 7nl n 2 
po \cj dz dz J 7 

If 7 > 1, then N§ > (subadiabatic thermal stratification) 
and the equilibrium vertical structure of the disc is convec- 
tively stable. For 7 < 1, then N 2 < (superadiabatic ther- 
mal stratification) corresponding to a convectively unstable 
equilibrium. For 7 = 1 (adiabatic thermal stratification), 
No = and all motions/modes due to buoyancy disappear. 
Since our goal here is to investigate the linear coupling be- 
tween vertical convection and SDWs and the possible role of 
this phenomenon in the disc dynamics and angular momen- 
tum transport, we focus, as in RG92, only on the superadi- 
abatic vertical structure choosing 7 = 0.8 throughout this 
paper. This value of 7 is close to unity reflecting the fact 
that the degree of superadiabati city is not very large in pro- 
toplanetary discs (|Rafikovll200"7l ). Strictly speaking, the case 
7 < 1 does not make much thermodynamic sense, because 
it implies c v > c p . In spite of this, there is no problem with 
equation (3) describing entropy conservation. In principle, 
the convectively unstable regime can be modelled without 
using the condition 7 < 1, if we appropriately choose cool- 
ing and heating functions in the entropy equation (see also 
RG92). So, a simple 7-prescription mimics the basic features 
of convectively stable/ unstable discs well without introduc- 
ing additional heating and cooling functions that might un- 
necessarily complicate the situation. 



2.2 Perturbation equations 

Consider now small perturbations to the equilibrium state 
(8) with the background flow uo. Linearising equations (1-3) 
about this state, we obtain the following system governing 
the perturbation dynamics 



Du' x 
Dt 



po ox 



Du 'v 1 dp' , 
Dt po ay 



Du' z _ 
Dt 

Dp' 



po dz 
V • (poll') 



9o- 



po 



Dt 

Ds' 7^0 , 

— — u, 

Dt g z 



(9) 
(10) 

(11) 
(12) 
(13) 
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where 

D d _ d 

and u', p , p' , s' are the perturbed velocity relative to the 
background Keplerian shear flow, perturbed density, pres- 
sure and entropy, respectively. From the definition of specific 
entropy above we also have the relation 



P_ 

Po 



7- 



Po 



(14) 



The form of equations (9-13) permits a decomposition of the 
perturbed quantities into shearing plane waves, or spatial 
Fourier harmonics (SFH) with time-dependent amplitudes 
and phases, 



A4(r,t)\ 






u' y (r,t) 




Uy(t) 


u' z (v,t) 




U z (t) 


\s'(v,t)J 




V *(*) / 



exp 



J£l_ 
2H 



P'(r,t) 
p'(r,t) 



exp 



2H 



exp(ik x (t)x+ik y y+ik z z), 
(15) 

exp (ik x (t) x + iky y + ik z z) , 



(16) 



k x (t) = qQkyt. 



changed, whereas the radial wavenumber k x (t) varies with 
time at a constant rate qQk y if k y 7^ (i.e., for non- 
axisymmetric perturbations) due to sweeping of wave crests 
by the background shear flow. For convenience, in k x (t) we 
have shifted the origin of time towards negative values, so 
that for t = 0, k x (0) = as well. In this case, an initially 
leading SFH with k x (t)/k y < at t < eventually be- 
comes trailing with k x (t)/k y > at t > 0. This change of 
SFH's orientation from leading to trailing is called 'swing' 
and occurs when k x (t) — 0. This technique of decomposition 
of perturbations into shearing wav es was originally devised 
by Lord Kelvin (|Thompsonl fl887) in order to study tran- 
siently growing solutions in inviscid incompressible parallel 
flows with linear shear. It is widely used today in various 
applications involving flows with a linear shear of velocity 
profile and greatly helps to grasp intermediate-time linear 
phenomena - transient growth and coupling of perturba- 
tion modes - in shear flows, which tend to be overlooked in 
the standard modal analysis. Besides, such shearing waves 
are really the only spectral basis compatible with shearing- 
periodic boundary conditions, which are almo st always used 
in loc al simulations of accretion discs (see e.g.. lllawlev et all 
1995). We would also like to mention that rec ently it has 
been mathematically proven by Yoshidal (|2005l ) that shear- 
ing waves in fact represent the simplest/basic 'elements' of 
dynamical processes at linear shear. The exponential fac- 
tors in (15) and (16) involving ±\z\/2H are necessary in 
order to compensate for otherwise exponentially increasing 
with height perturbation energy due to stratified equilib- 
rium. With the latter form, the perturbation energy, propor- 
tional to pou' 2 , is vertically uniform when averaged over the 
vertical wavelength (see also lLerche fc Parkerlll967l . RG92). 
The inclusion of these exponential factors also ensure that 
there are no imaginary terms in the dispersion relation de- 



rived below that could lead us to infer some sort of spurious 
instability. 

Substituting (15) and (16) into equations (9-14) and 
using replacement (7), we arrive at the following system of 
first order ordinary differential equations that govern the 
linear dynamics of SFHs of perturbations 

du x 



— - -ik x (t) — + 2Qu y , 
at pm 



dp 



du z _ 
~dT ~ 

= — pm 



diiy 
~dT 

ik z 



iky 

Pr 

sign(z) 
2H 



+ {q-2)Qu x , 



f ■ f \ f 

sign(z)p , 

Pm pm 



ik x {t)u x + ikyUy + I ik 



ds 
~di 



sign(z) 
2H 



-sign(z) u z 



(17) 
(18) 
(19) 
(20) 
(21) 



■ 7- 



(22) 



Pm Pm 

Note that by means of the exponential factors in (15) and 
(16), we have rendered the coefficients in these equations z- 
independent, though they still contain the sign of z. Below 
we will reduce this set so that sign-dependence will disap- 
pear too. 

It can be readily shown that the system (17-21) pos- 
sesses an important time invariant - the linearised potential 
vorticity perturbation 



/ = ik x (t)u y 
(2 - q)Q 



P 

pm 



sign(2p 
7iV 2 



sign(z) 
2H 



const, 



whose conservation is a direct consequence of Ertel's poten- 
tial vorticity theorem (4). This time invariant /, in turn, 
indicates the existence of the vortical/aperiodic mode in 
the perturbation spectrum, which is characterised by non- 
zero potential vorticity (/ 7^ 0) and represents a station- 
ary (i.e., time-independent) solution of equations (17-21) in 
the absence of shear. In the present three-dimensional case, 
this vortical mode originates from the combined action of 
the vertical gravity force (stratification) and the Coriolis 
force and represent s a 3D gene r alisat io n of the 2D vortical 
mode c o nsidered in Bodo et a l. (2005); Joh nson fc Gammiel 
i|2005aft ; lMamatsashvili fc Chagehshvilil l|2007l ); HP09a. Be- 
cause of the shear/non- normality of disc flow, the vorti- 
cal mode can undergo large transient amplification and 
effectively (linearly) couple with and excite other modes 
(SDWs, inertia-gr avity waves, bar oclinic modes, etc.) ex- 



isting in the disc ( Tcvzadzc ct al. 2003; Bodo et al.l 



Mamatsa shvili fc Ch agclishvili 2007; Tcvz adze et al.l 



2005; 



2008, 



201(1 HP09a). However, in the linear regime considered here 



with no disc self-gravity, the transient (algebraic) growth of 
the vortical mode is overwhelmed by the exponential am- 
plification of the convective mode, so we exclude the former 
mode from our analysis by setting the linearised potential 
vorticity to zero, / = 0. In the non-linear regime, however, 
the exponential growth of the convective mode is saturated, 
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so convection no longer dominates over vortices and both 
can equally take part in the SDW generation process. 

Taking into account that by our choice I — 0, we can 
eliminate u x , ii y , u z and p in favour of p and s in equations 
(17-22). Then, for convenience, switching to a new auxiliary 
quantity 



I = 



1 



1 sign(z)ff 
7 -yNg 



f>n 



ik z — 



sign(z) 



2H 



7^ 2 



sign(z) 
2H 



instead of the density p, we finally arrive at the second or- 
der system with real coefficients forming the basis for our 
subsequent analysis, 



d I 

= -[C 2 s k±(t) + K 2 (t)]l - cjk ± (t)s, 



dfs _ 
where 



-1.) 

4H 2 ) 



kl + 



(23) 
(24) 



K 2 (t) = 2(2 - q)tt 2 



4qVL 2 k 2 y 3q 2 Q 2 k* 

~kJW + ~kJW' 



and k±(t) = k 2 x (t) + k 2 = fc 2 (l + q 2 Q 2 t 2 ). Although the 
physical meaning of I is less transparent, by taking its ad- 
vantage we have reduced the original system (17-22) to a 
more compact one (23) and (24), with all the coefficients be- 
ing real, in order to ease further manipulations and to make 
our equations friendlier for numerical treatment. Equations 
(23) and (24) contain all the information on various pertur- 
bation modes (apart from the vortical mode) present in a 
vertically stratified Keplerian disc. In the next section, we 
remove shear from these equations for the moment in order 
to classify modes in the disc and then put shear back again 
to see how it alters the mode dynamics and what new effects 
it introduces. 



2.3 Classification of modes in the absence of shear 

Consider first a simple case without shear (i.e., a rigidly 
rotating disc) by setting q — in equations (23) and (24). 
After that, all the coefficients in these equations become 
time-independent and, therefore, we can look for solutions 
in the form I, s oc exp(iLut). Substituting this into equations 
(23) and (24), we obtain the following dispersion relation in 
the absence of shear 



2,2. 21,2, 1 \ . 

c s k ± +c s [ k z + I + k 



Nokl + k 2 lk 2 z + 



4H 2 J 



0, 



where k 2 — 4f2 2 is constant when q = 0. This dispersion 
relation has two different solutions corresponding to two dif- 
ferent types of perturbation modesQ 



If we included the vortical mode, it would correspond, as noted 
above, to the stationary solution u> = 0, since 1^0 for this mode. 



1. A high-frequency SDW mode with 
2 _ c 2 k 2 ± +c 2 (k 2 z + l/4H 2 ) + K 2 



l + Wl 



4c![7V 2 fci + K 2 (fcf + 1/4H 2 )] 
[clfc 2 + c|(fcl + l/4iP) + K 2 ] 2 



(25) 



the restoring force for which is mainly provided by com- 
pressibility/pressure forces, but is modified by stratifica- 
tion/bouyancy and the Coriolis force. So, in this mode, 
pressure perturbations dominate over entropy perturba- 
tions. For large wavenumbers k±H,k z H 3> 1 (i.e., for 
wavelengths much smaller than the disc scale height), the 
effects of stratification and rotation are negligible and the 
frequency of the SDW mode reduces to uj 2 ~ C^{k\ + k 2 ). 
In an unstratified disc, for z— independent perturbations 
(i.e., in the razor-thin disc approximation), from (25) we 
obtain lj 2 — c 2 (k 2 + k 2 ,) + k 2 , which is a classical dis- 
persion relation of two-di mensional SDW s in non-self- 
gravitating discs (see e.g., iBalbusI I2003T ) and therefore 
expression (25) can be regarded as a generalised form 
of the SDW mode frequency in the presence of vertical 
stratification. 
2. A low-frequency inertia-buoyancy mode with 



2 

W s = 



\ki +c 2 s (k 2 + 1/4H 2 ) + k 2 



1- 4/1 



4c 2 [N 2 k ,2 _ + n 2 (k 2 + 1/4H 2 )] 



\c 2 k\ + c 2 (k 2 + 1/4H 2 ) + k 2 



(26) 



the restoring force for which is mainly provided by 
vertical buoyancy and the Coriolis force, but modified 
by compressibility. So, in this mode, entropy perturba- 
tions dominate over pressure perturbations. The inertia- 
buoyancy mode represents an inertia-gravity wave in 
the case of subadiabatic and adiabatic stratifications 
(N 2 ^ 0) and the convective mode in the case of 
superadiabatic stratification (Nq < 0), which is con- 
sidered here(f| The dynamics of inertia-gravity waves 
in K eplerian discs is extens i vely studied elsewhere (see 
e.g.. iLubow fc Pringlelll993l : IOgilvielll998l: IBalbusI 1200 "' 



iTevzadze et alj 120031 . I2OO8I ; lLatter fc Balbud l2009h and 



will not be dealt with here. Again, in the limit of large 
wavenumbers k±H,k z H 2> 1, the effect of compressibil- 
ity on the inertia-buoyancy mode is small and we get the 



well-known dispersion relation cj 



which 



can also be deri ved using the incompr essible (Boussinesq) 
approximation l)Tevzadze et al.ll2008l ). 

Thus, in a stratified compressible disc, perturbations can 
be classified into two basic types - the SDW and inertia- 
buoyancy modes (apart from the vortical mode) and any 
general perturbation (with zero potential vorticity) can be 
decomposed into the sum of these two modes. We should 
note, however, that the mode classification performed here is 
actually applicable in the case of the constant-g approxima- 
tion (7), which is, strictly speaking, valid when characteristic 

2 Although we call this a convective mode, it actually exhibits 
convective instability (i.e., it has lo 2 < 0) only when the wavenum- 
bers satisfy K 2 (k 2 + 1/4H 2 ) < — A^q/c^ in order to overcome the 
stabilising effect of rotation. 
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vertical length-scales of perturbations are less than the disc 
scale height. A more general classification of vertical modes 
in stratified discs without invok i ng this approximation i s 
performed inlRuden etafl (I1988T); iLubow fc Pringld (|l993h : 
iKorvcanskv fc Pringld (|l995h ; lOgilvid (|l998l ) fin these pa- 
pers, the compressible p-mode corresponds to our SDW 
mode and the r- and convective g-modes to the inertia- 
buoyancy mode). Nevertheless, the classification adopted 
here allows us to grasp the key effects of shear on the dynam- 
ics of the inertia-buoyancy (convective) and SDW modes. 

To analyse the behaviour of these two modes, we define 
new characteristic quantities, or eigenfunctions, ip s and tp g , 
for them as 



■0s 



k ± [c 2 s (k 2 z + 1/4H 2 ) - N%}1 - [a; 2 - c 2 (fc 2 + l/W 2 



(27) 



, _ k ± [c 2 (k 2 z + 1/AH 2 ) - N$]l - [u£ - c 2 (k 2 z + l/AH 2 )]s 

ft — o 2 

- lu£ 

(28) 

describing, respectively, the SDW and inertia-buoyancy 
modes. These eigenfunctions are convenient and physically 
revealing, as substituting them into equations (23) and (24) 
with shear terms removed, we obtain two separate equations 
for each mode eigenfunction 



d 2 ip s . 2 , 



0, 



dt 2 



+ UJgl/jg = 



(29) 



that allow us to study these modes individually. All the 
other quantities (u x ,u y ,u z , p,p) can be expressed through 
the mode eigenfunctions il> a ,il>g and their respective first or- 
der time-derivatives. Hence, we can fully determine the per- 
turbation field corresponding to a specific mode by setting 
the eigenfunction of the other mode and its time deriva- 
tive to zero. In other words, if we want to have either only 
SDWs or the inertia-buoyancy mode in the disc flow, we 
should initially set simply tp g (— oo) = dip g /dt(—oo) = 
or Vs( — °°) = dip s /dt(—oo) = 0, respectively. Modal equa- 
tions (29) governing the dynamics of the SDW and inertia- 
buoyancy modes are not coupled, implying that in the ab- 
sence of shear the perturbation modes evolve independently, 
that is, initially exciting one either mode with a specific char- 
acteristic time-scale does not lead to the excitation of other 
modes with different time-scales. In the following sections, 
we investigate how shear modifies modal equations (29) 
and its influence on the dynamics of the SDW and inertia- 
buoyancy modes. Specifically, we will show that in the pres- 
ence of Keplerian shear, these equations are no longer inde- 
pendent from each other, i.e., become coupled that, in turn, 
results in the linear coupling between these two modes and, 
in particular, between SDWs and convection. 



and length-scales (wavenumbers) as well as I using c 3 /fi as 
a unit of length, 

(K x ,K y ,K z ) = ^-(k x ,k y ,k z ), — H = —^1 
il c s ^7 c„ 

(s,ip a> ipg are already non-dimensional). Accordingly, the 
normalised radial wavenumber K x varies with normalised 
time as K x {t) = qK y t. Also, as defined above, K±(t) — 
K x (t) + K 2 = K 2 (l + qt 2 ). Since K x (t) and K 2 (t) in equa- 
tions (25) and (26) are now time-dependent as a result of 
shear, the non-dimensional frequencies uj s and uj g also be- 
come functions of time: 

ull{t)= ^) + KUt) + Kl +1 ,^ 



i + Wi 



4[N 2 Kl(t) + ^(t)(KI +J /4)] 
[n 2 (t) + Kl{t) + K 2 + -y/4] 2 



K 2 (f) + ^(t) + if z 2 + 7 /4. 



A[N 2 Kl(t) + K 2 (t)(KI+j/^] 
[ K 2 (t) + K'i(t) + K 2 + J /4] 2 



At large times, K±(t) — > oo and from these expressions it 
follows that co 2 (t) ~ K\(t), to g (t) ~ Nq. As noted above, all 
the other quantities can readily be expressed through ip a ,ipg 
and their first order time derivatives, so we carry out the 
further analysis working mostly with the eigenfunctions. 

Expressing I, s through ips,tpg from equations (27) and 
(28) and substituting them into equations (23) and (24), 
we arrive at the following system of coupled second order 
differential equations for the eigenfunctions 



d 2 jj a 
dt 2 



d 2 % 



+ fs^ + (^ + Ac 2 )V> s = f 9 ^+ A^ 2 ft , (30) 
+ + K 2 + ^ 9 2 )ft = + (31) 



where the new coefficients f s , Aw 2 and f g , Auj 2 , compared 
with equations (29), describe modifications to the dynam- 
ics of the SDW and inertia-buoyancy modes and their cou- 
pling brought about by the shear of disc flow. They depend 
on time as well as on the azimuthal K y and vertical K z 
wavenumbers: 



f 3 {Ky,K Z ,t) 



f g (Ky,K Z ,t) 



2K ± (t) / Ws 2 (t)-JTf- 7 /4 



w|(t)-w|(t) 



K x {t) 



2K ± (t) fuj 2 (t)-K 2 -j/A 



w|(t)-w?(t)V K±(t) 



2.4 Effects of shear on mode dynamics 

Here we derive the generalised modal equations from equa- 
tions (23) and (24) including shear q. For the sake of further 
analysis, let us first non-dimensionalise the basic quantities 
by obvious units: time (frequencies) by means of the angular 
speed Q, 



fit 



K 



Q 2 



N 2 



Aoj 2 s (K y ,K z ,t) 



Auj 2 g (K y ,K z ,t) 



K ± [t) (u 2 {t)-K 2 z - 1 /A 



w2(t)-wg(t)V K x (t) 



K ± (t) fu 2 g (t)-K 2 - 7 /4 



wa(t)-w2(t)V K ± (t) 



(here and below primes denote the time derivative). Note 
that all the time derivatives in these coefficients are pro- 
portional to the shear parameter q as well as to the 
azimuthal wavenumber K y (time enters through qK y t) 
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and thus vanish in the shearless limit q = or in the 
case of axisymmetric perturbations with K y = 0, re- 
ducing equations (30) and (31) to decoupled equations 
(29). In vertically stratified Keplerian discs, the dynam- 
ics of axisymmetric pertu rbations (modes) was exten- 
sively studied in the past dRuden et al.l [l988l; [Kiev et all 



sively studied m tne p ast [rtudc n et al.l lli JBg; rucy et ai, 
1993 : iLubow fc Pringld ll993TlKorvcanskv fc Pringlel Tl995l : 



Qgilvid Il998l ; iMamatsashvili fc Ricel l201of T~ so we omit 
them from the present analysis and concentrate only on 
non-axisymmetric perturbations having non-zero azimuthal 
wavenumber K y 7^ 0, which we assume to be positive 
{K y > 0) throughout without loss of generality. 

Equations (30) and (31) describe the linear dynamics 
of the SDW and inertia-buoyancy modes and their coupling 
in a vertically stratified, compressible disc with (Keplerian) 
differential rotation. As mentioned above, novel features in 
the mode dynamics, in comparison with the case of no shear, 
arise from the coefficients f s , Auj 2 and f g , Auj 2 that originate 
from shear. The homogeneous (left-hand side) parts of equa- 
tions (30) and (31) describe the individual behaviour of the 
SDW and inertia-buoyancy modes, respectively, in the pres- 
ence of shear. The most interesting aspect here is the source 
terms on the right-hand side of these equations that provide 
coupling between these two modes. In other words, initially 
imposed only one mode acts as a source for another mode 
and can excite it in the course of evolution. Thus, the shear 
of disc flow introduces a phenomenon of mode conversion. 
We will see below that because of such a coupling, the con- 
vective and SDW modes actually have strictly separate iden- 
tities only at large times (\t\ 3> 1), when the characteristic 
time-scales (frequencies/growth rates) of these modes dif- 
fer significantly (Fig. 1). As seen from Fig. 1, for finite times 
(\t\ £ 1) instead the mode time-scales are comparable and, as 
a result, the modes cannot be quite disentangled from each 
other, so that we have some mixture - 'convective-SDW - 
mode at such times. If there was not such a difference in the 
SDW and convective modes' time-scales, strictly speaking, it 
also would not make much sense to talk about homogeneous 
and particular inhomogeneous solutions of equations (30) 
and (31) separately, which themselves constitute a homoge- 
neous system. Below, we will demonstrate using numerical 
analysis how the linear mode coupling occurs in practice. 

For the purpose of numerical integration presented be- 
low, we change the eigenfunctions to tys = ip s /h s , *l/ 9 = 
ipg/hg, where 



h. 



cxp [—- f fsdt' j , h g = exp 



fgdt' 



These new eigenfunctions are more convenient, because sub- 
stituting them into equations (30) and (31), we arrive at a 
simpler system without first derivatives on the left hand side 

^ r +uu 2 ^ s = x S9 i^ + Xs 9 2* 9 , (32) 



d 2 t> 

~~d¥ 



Equations (32) and (33) with time-dependent (modified) fre- 
quencies 

- 2 2 . A 2 fs fs 
LO, = U}, + Acj„ — , 




q 2 = 2.25 



-0.5 
t 



Figure 1. Modified frequencies Co 2 and aj 2 as functions of time at 
K y = 2, K z = 3 for the convectively unstable stratification with 
7 = 0.8. For comparison, with dashed line we also show the square 
of the inverse shear (non-dimensional) time q 2 = 2.25. In the 
vicinity of t = 0, all these three time-scales become comparable, 
whereas at large times they get separated: ui 2 increases as oc t 2 
and ui 2 tends to constant Nq = —0.25. 



LU n = LU + AOJ - — f, 



- 2 
J 9 



Ik 

2 



respectively, for the SDW and inertia-buoyancy modes, re- 
semble equations of two coupled oscillators. These frequen- 
cies differ from corresponding cj 2 and w 2 due to the presence 
of the shear-induced terms. For large times as well as for 
K y ^> 1 and/or K z ^> 1, these terms are small and ii 2 and 
ujg go to corresponding values in the shearless limit: 

& a .{t) ~ Kl{t) + K$, 

The x parameters, describing the coupling between the 
modes, are given by 



Xsgl 



fghg 



f a h 3 



Xsg2 



ha I A,- 2 



f 2 
Jg 



— \ Aw., - — 

h a 



The coupling parameters Xsgi,Xsg2 describe the excita- 
tion of the SDW mode by the inertia-buoyancy mode, while 
Xgsi,Xgs2 describe the excitation of the inertia-buoyancy 
mode by the SDW mode. Figure 2 shows the temporal vari- 
ation of the coupling parameters during the swing of a per- 
turbation SFH from leading to trailing in the case of su- 
peradiabatic stratification with our chosen 7 = 0.8. Accord- 
ingly, from now on we speak only of the convective mode 
instead of the more general inertia-buoyancy mode. In this 
case, &g(t) has a negative sign for most of the time (Fig. 
1) and determines the shear-modified instantaneous growth 
rate of convective instability. It is seen from Fig. 2 that for 
fixed K y and K z , these parameters reach their maximal val- 
ues in the interval \t\ ^ 1, when the time-dependent radial 
wavenumber is not large, \K x (t) / K y \ = q\t\ 1, and rapidly 
decay at \t\ 3> 1, when \K x (t) / K y \ — q\t\ 1 as well. Thus, 
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Figure 2. Coupling parameters \ as a function of time for K y = 
2, K z = 3. They reach the highest values during \t\ H 1 and fall 
off at large times \t\ 3> 1. 
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Figure 3. Maximum values of modulus of coupling parameters, 
|x|j as a function of K y and K z . max|xs 9 i| and max|xs 9 2| are 
appreciable over a broader range of wavenumbers than max |x 9 sl| 
and max|x 9 s2|- 



because the coupling parameters \ evidently also originate 
from shear, we can conclude that at given K y and K z , the 
influence of shear on the mode dynamics can be important 
only for moderate radial wavenumbers and an efficient en- 
ergy exchange between the modes and the mean disc flow 
should be expected to occur just during this interval, which 
we will call a coupling interval. For \K x (t) / K y \ 3> 1, the 
coupling parameters are small and hence the modes become 
dynamically decoupled from each other and evolve indepen- 
dently (see the next section). In Fig. 3, the maximal values 
(over the interval \K x (t) / K y \ 1) of the modulus of the 
coupling parameters are plotted in the (K y , K z )- plane. The 
maximum values of |xs 9 2|, |Xssi| and |x 9 s2| are appreciable 
for of the order of unity or smaller K z , K y and rapidly decay 
with the increase of these wavenumbers. The maximum of 



|Xs 9 i| is appreciable over a broader range of wavenumbers, 
but it has larger values at smaller K y . Because the coupling 
parameters are generally one of the factors determining the 
generation of one mode by another, we should expect the 
efficiency of coupling of the SDW and convective modes as 
a function of K y and K z to follow, to a certain extent, the 
same trend as that of the coupling parameters shown in 
Fig. 3. Although in a superadiabatically stratified disc the 
generation (triggering) of convective motions by SDWs, de- 
termined by the coupling parameters x 9 siiX 9 s2, is in prin- 
ciple possible, here we focus only on excitation of SDWs by 
convection, determined by Xagi, Xsg2, as it seems more inter- 
esting. Besides, convection can be easily triggered by other 
mechanisms in discs. 



3 SDW AND CONVECTIVE MODES - 
SHEAR-INDUCED COUPLING 

In this section, we study the shear-induced dynamics of the 
non-axisymmetric SDW and convective modes by numeri- 
cally solving equations (32) and (33). At large times/radial 
wavenumbers (\K x (t)/K y \ — q\t\ 2> 1), the adiabatic 
(WKBJ) condition with respect to time is satisfied for the 
mode frequency /growth rate, 

I « l<tf, 9 (t)l- 

A physical interpretation of this inequality is that the time- 
scale of shear-induced variation of the frequency /growth rate 
is much larger than the characteristic time-scale (i.e., inverse 
frequency /growth rate) of the mode itself. In this adiabatic 
regime, o)g(t) ~ Nq < 0, so it is nearly independent of time 
and is much less than the frequency of SDWs, which in- 
creases linearly with time, u) 3 (t) ~ K(t) sa qK y \t\ 2> 1, so 
that the mode time-scales are well separated (Fig. 1). In ad- 
dition, all four coupling parameters x are negligible at these 
times (Fig. 2). These imply that in the adiabatic regime, 
shear plays only a minor role in the dynamics of SDWs and 
convective instability. Below we will refer to this asymptotic 
stage as adiabatic/non-coupling region in the K-space. Con- 
sequently, if the convective mode is in the adiabatic region, 
it will evolve with time as in the shearless limit, without any 
exchange of energy and exciting other modes (i.e., SDWs) 
in the disc. However, if K x (t) crosses the coupling interval 
| K x (t) /K y I ^ 1 during its drift, where, as we have seen, the 
influence of shear on the mode dynamics becomes signifi- 
cant, the excitation of SDWs by non-axisymmetric convec- 
tive motions can take place. So, the goal of this section is to 
examine this phenomenon in detail. 

For further use, we also introduce the mode energies as 

E s EE |*1| 2 +£ s 2 (i)|* s | 2 , E g = |* 9 | 2 +LJ 2 g (t)\* g \ 2 . 

for the SDW and convective modes, respectively. In the ab- 
sence of shear they are conserved quantities, as it follows 
from equations (29). The mode energies are useful in char- 
acterising mode dynamics for asymptotically large times, in 
the adiabatic regime. In this regime, the energies vary with 
time due solely to the mean shear flow. In particular, the 
SDW mode energy grows linearly with time, 

E s oc w 8 (t) ~ K(t) oc \t\, 

which follows from equation (32) in the WKBJ regime ne- 
glecting the coupling parameters, which are small in this 
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regime anyway. Such an asymptotic behaviour of the SDW 
mode energy will later prove to be a useful diagnostic in 
analysing the generation of SDWs. Note also that at large 
times E„ is obviously positive definite, whereas E g can be- 



come negative (Fig. 4), since tb g ~ No < at these times. 



3.1 Generation of the SDW mode by the 
convective mode 

We wish to investigate the convective instability in the pres- 
ence of Keplerian shear and, specifically, how convective mo- 
tions can excite SDWs. To this end, initially at t — —to, 
where to 2> 1 is some large positive parameter, we impose 
a tightly leading (i.e., with K x (-t )/K y < -1) SFH of 
the convective mode on the flow without any mixture of 
a SDW mode SFH and follow the subsequent evolution of 
the eigenfunctions and perturbed quantities until t — to, 
at which stage their SFHs become tightly trailing (with 
K x (to) / K y 3> 1). For numerical integration we use a stan- 
dard Runge-Kutta scheme (MATLAB ode45 RK implemen- 
tation). 

To prepare such initial conditions, we make use of the 
fact that at t <C —1 the dynamics is adiabatic and the modes 
do not interact with each other that, in turn, permits us to 
pick out (impose) only the convective mode at the begin- 
ning. In this adiabatic/non-coupling regime, the convective 
mode is given by the WKBJ asymptotic solution of the ho- 
mogeneous part of equation (33), which we take to have the 
form 



* 



Co 



■exp 



u g (t')\dt' ) at t < -1, (34) 



where Co is a some arbitrary constant setting the initial 
value of *i!g at the start of integration at t = —to- This so- 
lution means that the convective mode is evanescent in the 
distant past (at t — > — oo) and grows with time afterwards. 
As for the initial value of the SDW mode eigenfunction, a 
high-frequency SDW component is absent at the outset, so 
that ty s should be given only by non-oscillatory particular 
solution, \&i s ', of equation (32) coming from the source term 
on the right hand side that is associated with the convective 
mode. Asymptotically at t <C — 1, we can ignore the second 
order time derivative in equation (32) and represent the ini- 
tial value of ty s at these times, to a good approximation, 
as 



(35) 



Indeed, as mentioned above, in the adiabatic regime, (i) 
the period of the SDW oscillations is much smaller than 
the shear/dynamical time, or equivalently \£j' s \ <C w,, (ii) 
the coupling parameters vary a little during the oscilla- 
tion period and at the same time (iii) the frequency of 
SDWs is much larger than the characteristic growth rate 
of the convective mode, \ui g \ < ii s . As a result, the sec- 
ond time derivative of the particular solution (35) turns out 
to be much smaller than the right hand side term of equa- 
tion (32), thereby validating our approximation. This slowly 
varying solution ensures that the oscillatory SDW mode is 
absent in the initial conditions. Note also that as the cou- 
pling coefficients are vanishing at t — > — oo, we also have 
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Figure 4. Evolution of the eigenfunctions *S> g and entropy s, 
normalised perturbed pressure p/pm and corresponding energies 
Eg , E s pertaining to an initially imposed purely convective mode 
SFH with K y = 2, K z = 3. At the beginning, ty g and s evolve adi- 
abatically, exponentially growing with the characteristic growth 
rate \uig\, while ty s and p are still nearly zero. At around t = 0, 
rapid oscillations abruptly emerge in the evolution of ^ s and p 
that indicate a trailing SFH of the SDW mode being excited. In 
addition, the convective mode causes p to grow as well while os- 
cillating. Accordingly, the energy of the SDW mode, E s , being 
negligible at t < 0, after the transient amplification event in the 
vicinity of t = 0, increases linearly with time in the subsequent 
adiabatic regime at t 2> 1 (which actually starts already at t ~ 1). 



the condition of the absence of SDWs in the shearless limit. 
Thus, asymptotic solutions (34) and (35) can be viewed as 
a generalisation of the recipe for imposing only the convec- 
tive mode, without a mixture of SDWs, on the flow in the 
absence of shear (i.e., 1 3/ s (— oo) = ^f' a (— oo) = 0) to the case 
of non-zero shear provided, however, that the adiabatic ap- 
proximation is still met, which is always the case at large 
times. 

Figure 4 shows the subsequent temporal evolution cor- 
responding to initial conditions /solutions (34) and (35). To- 
gether with the eigenfunctions, which more clearly illustrate 
mode coupling, we also plot the time-development of the per- 
turbed entropy s, normalised pressure p/p m and mode ener- 
gies Eg and E a . At the beginning, while the convective mode 
SFH is still in the adiabatic region K x (t) /K y — qt <C — 1 and 
far from the coupling region, the eigenfunction ^ g , following 
expression (34), grows exponentially with the instantaneous 
growth rate \ui g \ and induces a similar behaviour in the en- 
tropy, whereas evolving according to (35), still remains 
very small and non-oscillatory. The evolution of the corre- 
sponding energy E„ is similar to that of >I/ S . It is also seen 
in Fig. 4 that at this initial adiabatic stage of evolution, the 
entropy perturbation is dominant over the pressure pertur- 
bation, which, like ^s, is small and aperiodic initially. It 
can also be shown that the form of our initial conditions for 
the convective mode implies E g = at t <C — 1, because of 
the negative u)g < 0; of course, this does not mean that the 
convective mode is zero. 

As time passes, the time-dependent radial wavenumber 
K x (t) of the convective mode SFH gradually approaches the 
coupling region \K x (t) / K y \ 1 and begins to cross it. In 
this region, the effects of shear come into play: the coupling 
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parameters become appreciable (Fig. 2) and simultaneously 
the characteristic time-scales of the SDW mode, convective 
mode and the shear time (= 1/q in the non-dimensional 
form) become comparable |f| or equivalently q ~ |tD s | ~ ui s 
(Fig. 1), also implying that in the coupling region, the 
adiabatic condition for the SDW frequency breaks down, 
\u)' s \ ~ \u>g\, that is, the coupling region is also non-adiabatic. 
This results in an efficient energy exchange between these 
two modes and the background shear flow in this coupling 
region. First, at — f t < 0, the eigenfunction *I/ S , pressure 
p and the energies E g ,E s undergo transient amplification 
by extracting energy mainly from the background flow, but 
still remain aperiodic. The evolution of ty s has now devi- 
ated from the asymptotic solution (35), because the second 
derivative of \t s has become important in equation (32), 
though it is still following non-oscillatory vfi 8 . Then, during 
a short period of time around t — 0, when K x (t) crosses the 
point K x — 0, swinging from negative (leading) to positive 
(trailing), rapid oscillations abruptly appear in the evolu- 
tion of ty s and the pressure, signaling the generation of a 
trailing SFH of the SDW mode. So, now at t > 0, there are 
the newly excited trailing SDW mode SFH and the former 
convective mode SFH. It can be said that the convective 
mode, in some sense, acts as a mediator between SDWs and 
the disc shear flow. The energy needed for the SDW exci- 
tation is mainly extracted from the shear with the help of 
the convective mode. We would also like to note that, as 
mentioned above, because the time-scales of the convective 
and excited SDW modes are comparable, strictly speaking, 
these modes cannot be well distinguished from each other 
in the coupling region. In other words, we have a mixture 
- 'convective-SDW mode - at intermediate times, \t\ ^ 1. 
Then, on leaving the coupling region, K x (t) moves into the 
next adiabatic region K x (t)/K y 3> 1 where the linear dy- 
namics of SDWs and convection become decoupled with no 
further energy exchange, so the modes thus acquire truly 
separate identities. The mode time-scales has now been sep- 
arated (Fig. f): the pressure p rapidly oscillates as a result 
of the new SDW component in $ s with the frequency &s(t), 
which increases linearly with time and soon becomes much 
larger than the growth rate, |i£ fl j, of the convective instabil- 
ity. The eigenfunction tyg continues to grow exponentially 
with this growth rate and causes a similar growth in the 
entropy. As typical of SDWs in the adiabatic regime, the 
corresponding energy E 3 increases linearly with time solely 
due to the background flow. 

Finally, we would like to stress again that the stud- 
ied here linear coupling of the SDW and convective modes 
caused by shear is essentially of the same nature as other lin- 
ear mode coupling phenomena occurring in disc flows with 
Keplerian shear: coupling of vortical and wave modes, that 
is, the generation o f SDW s by vortices dBodo et all 120051; 
Johnson fc Gammid l2005al ; iMamatsashvili fc Chagelishvilil 
20071. HP09a), coupling of v ortices and inertia-gravity waves 
i Tevzadze et al_. 2003, 2008) and coupling of baroclinic and 
SDW modes (|Tevzadze et alj|20ld ). Thus, the linear mode 



3 This is not the case at large K y and/or K z , because the mode 
time-scales still remain separated even for \K x (t) / K y \ ^ 1 and 
therefore the mode coupling is negligible at large wavenumbers; 
see section 3.2. 



coupling is a generic phenomenon inevitably taking place 
whenever the mean flow is inhomogeneous (i.e., the veloc- 
ity profile has a non-zero shear) and all these cases are its 
special manifestations for particular physical settings. In 
this connection, we should mention that it also occurs in 
MHD shear flows, where in general there are larger number 
of mode branches with different time-s cales and coupling 
amon g them becom es quite complex ( Chagclishvili ^et al.l 
1996). For example, iHeinemann fc Papaloizoul (|2009bl ) de- 
scribed the generation of SDWs by the vortical mode in a 
MRI-turbulent disc shear flow. The coupling among differ- 
ent MHD modes can also be seen in the line a r ana lysis of 
non-axisymmetric MRI by iBalbus fc Hawlevl fl992), how- 
ever, the authors do not identify it as such. Thus, the pri- 
mary effect of shear of the disc's differential rotation on the 
perturbations dynamics is that it couples, or introduces a 
new channel of energy exchange among different types of 
perturbation modes existing in the disc as well as with the 
disc flow. In the present case, the linear dynamics of the 
convective mode is accompanied by the generation of high- 
frequency SDWs during a finite-time interval as SFHs swing 
from leading to trailing. As noted before, the shearing wave 
(non-modal) approach, because of not involving spectral ex- 
pansion in time, allows us to trace an entire temporal evo- 
lution of perturbations and thus to reveal new aspects of 
non-axisymmetric convection in disc shear flows. This, in 
turn, can provide further insight into the role of convection 
and its coupling with SDWs in angular momentum trans- 
port, especially now when the capability of convection to 
transport angular momentum outwards in discs has recently 
been established via numerical simulations (we address the 
transport properties of the modes in section 4). Previous lin- 
ear studies of non-axisymmetric convection have been either 
in the framework of the modal approach (|Lin et al.lll99St ) 
or using the shea r ing wave (non-modal) approac h (RG92; 
iKorvcanskvl 1 1993 : iBrandenburg k, Dintrana I2006T 1 as here, 
but without identification and characterisation of the mode 
coupling process. As pointed out in the Introduction, the 
results of the modal approach are actually applicable for 
asymptotically large times and therefore tend to overlook 
dynamical effects at intermediate times arising as a result of 
the flow non-normality/shear. Since the excitation of SDWs 
by convection is just one of such examples occurring during 
a limited time interval, in fact it cannot be captured in the 
framework of the modal approach. 



3.2 Generation of SDWs for various K y and K z 

Here we examine how the efficiency of the above-described 
SDW generation process depends on the azimuthal, K y , and 
vertical, K z , wavenumbers. Figure 5 shows the evolution of 
under the same initial conditions (34) and (35) as above, 
that is, consisting of an initially imposed only tightly leading 
convective mode SFH at various K y and K z , without a mix 
of the SDW mode SFH. For large K y and/or K z (i.e., when 
the azimuthal and/or vertical wavelengths of the modes are 
much smaller than the disc scale height), starting out with 
very small values, ^ s mainly undergoes only transient vari- 
ation in the coupling region (at \t\ ^ 1) followed by very 
weak, almost no generation of high-frequency oscillations, 
i.e., SDW component. The reason for such a behaviour is 
the following. At large K y and/or K z , the time-scales of 
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Figure 5. Evolution of ^ s corresponding to an initially imposed convective mode SFHs with various K y and K z . Appreciable oscillations 
in the time-development of ^ s appear and, therefore an efficient SDW mode excitation occurs, at K y , K z ^ 1. By contrast, at large K y 
and/or K z , * s undergoes only small transient variation (growth) in the interval \t\ ,$ 1 with very weak, almost no SDW generation. 



the SDW and convective modes remain well separated dur- 
ing an entire course of evolution, not only at large times, 
and the adiabatic condition for the SDW frequency holds. 
As a result, non-oscillatory particular solution (35), valid 
for the initial adiabatic stage of evolution at t <ti — 1, in 
fact continues to hold at all times to a good approximation 
(because its second time derivative always remains negligi- 
ble in equation 32) and, consequently, an efficient coupling 
between the SDW and convective modes is not feasible at 
large K y and/or K z . So, in Fig. 5, the time-development of 
\l/ s at these wavenumbers is given by this particular solu- 
tion which does not involve any oscillations. (To be more 
precise, the SDW mode is still being generated, but with a 
very small, as evident from the panels with K y = 5, K z = 
5; K y = 10, K z = 5; K y = VQ,K Z = 1; K y = 10, = 0.6 
of Fig. 5 and also from Fig. 6.) Thus, the regime of large 
azimuthal and/or vertical wavenumbers is a weak coupling 
regime, where shear plays only a minor role in the mode dy- 
namics. By contrast, for smaller K y ,K z 1 (i.e., for wave- 
lengths comparable to the disc scale height) the favourable 
conditions for mode coupling - violation of the adiabatic 
condition and comparable time-scales of the SDW and con- 
vective modes and the shear time - can occur in the vicinity 
of t = 0, loading to an efficient generation of the SDW mode 
by the convective one, as described in section 3.1. 

We can quantitatively characterise the mode coupling 
as follows. At t > 0, the total solution for \T/ S consists of the 



two components 

where ^i s ', as before, is the non-oscillatory, slowly varying 
particular solution of equation (32) due to the convective 
mode and *s W ' ) is an oscillatory component related to the 
SDW mode generated by this particular solution. So, the rel- 
ative intensity of the wave generation process can be quan- 
tified by comparing the value of ^i w ' to that of ^i s \ To do 
this, we introduce the parameter 

c _ max t>0 jVI>< w) (t)l 
£ ~ maxi >0 |*i g) (t)| ' 

which is the ratio of the maximal values of |^i w ' | and \^i s ^ | 
over time (typically, both these functions achieve their max- 
imal values in the vicinity of t = 0, as seen in Fig. 5). Ob- 
viously, with our initial conditions (34) and (35), this ratio 
does not depend on the arbitrary parameter Co and is a func- 
tion of K y and K z only, which is plotted in Fig. 6. It is seen 
from this figure that at large K y and/or K z , as expected, 
e £ 0.1 is small, implying that the non-oscillatory particular 
solution dominates over the oscillatory SDW mode solution, 
i.e., wave generation is weak. For K y ,K z £ 1, e ~ 1 im- 
plying that the wave component is comparable to the non- 
oscillatory solution and therefore the SDW generation is ap- 
preciable. This also confirms the situation in Fig. 5. 
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Figure 6. Ratio, e, of the maximum values of the SDW mode so- 
lution to its generator non-oscillatory particular solution of equa- 
tion (32) due to the convective mode, plotted as a function of K y 
and K z . 
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Figure 7. Absolute value of the amplitude, \A\, of SDWs gen- 
erated by the convective mode as a function of K y and K z . It 
achieves the largest value |A| maa; = 0.082 at K y , m = 1.2, K z ,m = 
1.3 and decreases at smaller and larger K y , K z , implying also the 
decrease in the mode coupling efficiency at these wavenumbers. 



3.3 Amplitudes of generated SDWs 

We can go further and calculate the amplitudes of SDWs 
generated by the convective mode, which also represent a 
measure of the coupling efficiency between these two modes. 
In this section, we choose the constant Co in the WKBJ so- 
lution (34), describing an initially tightly leading convective 
mode SFH at t <C — 1, such that to have 



* 



Co 



exp 



\u g (t')\dt' 



An advantage of this form is that the calculated below SDW 
amplitudes depend only on K y and K z and are independent 
of the choice of the starting point —to (as long as to 2> 
1). In this initial adiabatic regime, as before, fy s = , 
because SDWs are absent at the outset. Then, in the non- 
adiabatic region at |t| 1, the mode coupling is at work 
and subsequently, at around t = 0, the SDW mode SFH 
abruptly emerges. After crossing the non-adiabatic/coupling 
region, in the next adiabatic region at t 3> 1, there are the 
former convective mode SFH and the generated by it SDW 
mode SFH, both with tightly trailing orientation. So, now, 
as noted above, the full solution for ty 3 is the sum of the 
following parts: 

* s = *i w) + *i g) = 
A 



— _exp -i / <2> s (t )dt } + 

^/LJ s {t) \ Jo 



+ 



A* 



exp ^i j u) s (t')dt'J + 

+ X'' 1 *' + *'''*« at t»l, 



where the first two terms are oscillatory WKBJ solutions 
of the homogeneous part of equation (32) that correspond 
to the excited SDWs. Since the initial conditions are real, 
these solutions come in complex conjugate pairs with differ- 
ent signs of frequency and with amplitudes A and its com- 
plex conjugate A* . Thus, the SDW mode SFH generated by 
the convective mode is actually a superposition of two SFHs 



corresponding to SDWs propagating in opposite directions. 
In other words, convection always pairwise excites SDWs 
that propagate oppositely, similar to SDWs generated by 
the vortical mode in 2D discs (see e.g., HP09a). The third 
term is again the non-oscillatory slowly varying particular 
solution due to the convective mode. 

Figure 7 shows the dependence of the absolute value 
of the amplitude, \A\, on K y and K z . When either of these 
wavenumbers is large, \A\ is small, because, as shown above, 
the SDW generation/mode coupling is inefficient. With de- 
creasing K y and K z , as expected, the amplitude increases 
and attains its maximum values at K y , K z 1, because an 
appreciable SDW generation takes place at such wavenum- 
bers, as seen in Figs. 5 and 6. Note that K y , K z 1 is the 
region where the coupling parameters Xsgi,Xsg2 are appre- 
ciable too (see Fig. 3). Comparing Figs. 6 and 7, we also 
see that the amplitude \A\ and e, which is a measure of the 
mode coupling efficiency, behave with K y and K z more or 
less similarly. 



4 ANGULAR MOMENTUM TRANSPORT BY 
THE CONVECTIVE AND SDW MODES 

In this section, we analyse the angular momentum transport 
properties of the SDW and convective modes. Specifically, 
we are interested in the radial component of angular mo- 
mentum flux carried by these two modes, because it deter- 
mines the radial transport of mass and hence accretion onto 
the central star. The angular momentum conservation law 
for linear perturbations in the shearing box model follows 
from the invariance of the system under translations in the 
azimuthal y-direction. It can be derived by application of 
Noether's the orem to the second o rder Lagrangian density 
(RG2; HP09a; iNaravan et al.lll987l ) 



L : 



Po 
2 



£1 

Dt 



Po 
2 



(7-l)(V-£) 2 + 



dxj dxi 



+ 



— + qpoil £ XJ 
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in this direction and yields 



D 



c)L 



d£,i 



Dt \d(D£i/Dt) dy J dx 3 \d(d^/dxj) dy 



dL 



LS 



J 2 



0. 



where £ is the displacement vector related to the perturbed 
velocity by 



Dj x 
Dt 



Dt 



+ qil£ x , Uz = ~[)t 



and Sij is the Kronecker delta and the summation is as- 
sumed over the repeated indices i,j — 1,2,3 corresponding 
to (xi, X2, X3) = (x,y,z). (It can also be readily shown that 
this Lagrangian density through variation gives the origi- 
nal set of linear perturbation equations 9-13). Since we are 
mainly concerned with the radial ^-component of the angu- 
lar momentum flux, we can average this equation over both 
y- and z-coordinates. After averaging, the canonical angular 
momentu m density in brackets in the first term, as demon- 
strated bv lNaravan et al.l (|l987T ). actually coincides with the 
density of the radial component of true physical angular mo- 
mentum, but taken with the minus sign and without fiducial 
radius ro that does not play any role in the local analysis. So, 
we write for the radial component, B, of angular momentum 
density 



B 



I dL g& 
\d{DZi/Dt) dy 



where the angle brackets denote averaging over y and z. On 
similarly averaging the second term and taking it with the 
minus sign, we obtain the desired radial component of the 
angular momentum flux 



F x = 



/ dL 
\d(dii/da 



dy 



= -( p 



dy 



where the pressure perturbation p' (as used in equations 9- 
13, the prime is not time-derivative in this case) is related 
to the displacement vector through p — — poc 2 s V ■ £ + pog£ z - 
After that the angular momentum conservation takes more 
compact form 



dB 
~dt 



dF x 
dx 



0. 



For the angular momentum flux associated with an individ- 
ual SFH given by (15) and (16) after spatial averaging we 
get 

F x = ^K v ]p(t)t x (t)-£ x (t)p*(t)], 

where again p(t) and £ x (t) are the amplitudes of SFHs de- 
pending only on time with p*(t) and £ x (t) being their re- 
spective complex conjugates. The velocity and displacement 
amplitudes are related by u x (t) = d£ x (t)/dt. We normalise 
displacement £ as other scale- lengths above, Q£/c s — s> £, and 
pressure and flux as p/p m — > p, F x /p m — > F x . Since, we have 
two modes - SDWs and convection, we decompose p and £ x 
into two parts corresponding to these modes 

P = Ps+Pg, ix = ix,s + ix.g, 

where p 3 and £ XlS are related to the SDW component and, 
therefore, are rapidly oscillating in time with the frequency 
of the latter, while p g and £ x , g are related to the convective 
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Figure 8. Radial angular momentum fluxes associated with the 
convective mode, F x<g , and with the generated by it the SDW 
mode, F x s , for K y m = X.2,K zrn = 1.3. The jump of F x s and 
F x ,g in the immediate vicinity of t = is due to the abrupt 
emergence of the SDW mode. 



mode and vary on its corresponding time-scale, relatively 
slowly compared with the SDW mode. So, the total flux 

Fx = ^K v [{p s +&)(£,, + £,„)- 

- (£ x , s +£x,g)(pt +P*g)] 

is also of oscillatory type because of the presence of the SDW 
mode, but we can conveniently smooth it by averaging over 
the oscillation period much smaller than the growth time 
of the convective mode (we denote the time-averaging with 
angle brackets without subscript). As a result, we have 

F X = ^Ky(p 3 £ XiS - £ X ,sP* s ) + jKy{p g C ig - £ X ,gP*g). (36) 

Thus, all the cross-products vanish after time-averaging and 
the total angular momentum flux appears as a sum of an- 
gular momentum fluxes 



F X , S = ~Ky (ps£ X , S - £x,sP* s 



£, X ,gP g 



F X ,g — ^Ky(j)g^ x g - 



related, respectively, to the SDW and convective modes (we 
have omitted the angle brackets in F x>g , as it does not 
change much during the oscillation period). In an analo- 
gous calculation of the angular momentum transport by non- 
axisymmetric shearing waves, RG92 obtained a rapidly os- 
cillating, on average negative, angular momentum flux and 
attributed it only to the convective mode without separat- 
ing/analysing the contribution from SDWs. In fact, this con- 
tribution is present, because, as demonstrated here, SDWs 
are inevitably generated by the convective mode due to 
background shear. So, decomposition (36) allows us to anal- 
yse the angular momentum fluxes carried by each mode in- 
dividually and contrast them. 

In Fig. 8, we show the time-history of the fluxes F XtS 
and F XtS when a tightly leading purely convective mode SFH 
is imposed initially. At t < 0, before SDW mode generation, 
angular momentum transport is solely due to the convective 
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mode SFH and is still small and positive. Then, in the vicin- 
ity of t = 0, trailing SDW mode SFH emerges, giving rise to 
its associated flux F x ,a at t > 0, which initially increases and 
soon settles to a positive and constant value. Thus, trailing 
SDWs always transport angular momentum outwards. (A 
similar result was also obtained by HP09a for SDWs gener- 
ated by the vortical mode.) As for the angular momentum 
flux of the convective mode, during a short period of time 
after t = 0, it is still positive, but quickly switches to in- 
creasingly negative due to the exponential growth of $ s , 
and dominates over the positive flux of the SDW mode with 
time. Because of the relatively short time-scale of the wave 
generation process in the immediate vicinity of t = 0, the 
variations of F Xt3 and F x<9 als o appear sharp around t — 0. 
However, as demon strated by iLesur fc Qgilviel i|2010h and 
iKapylaet alll|20irj ). this flux of angular momentum associ- 
ated with the convective mode, being negative in the linear 
regime, can actually change sign (direction) in the non-linear 
regime. We also see from Fig. 8 that shortly after the SDW 
mode generation, at K x (t)H /2ir ;$ 1 (i.e., when the radial 
wavelength is still larger than or equal to the scale height), 
its corresponding flux, F x<3 , achieves a maximum value com- 
parable to, though still less by absolute value, than that 
of the convective mode, F x , g . As iHeinemann fc Papaloizoul 
(2009b) showed, just these values of the angular momentum 
flux of SDWs during the time when their time-dependent 
K x (t) are still in the range K x (t)H /2tt ^ 1, are important 
and determine transport in the non-linear regime, because 
SDWs with radial wavenumbers larger than this are subse- 
quently damped via shock formation in the trailing phaseQ 
Thus, SDWs generated by convection are, in principle, capa- 
ble of fully contributing to the transport before undergoing 
significant damping. Based on this, in the non-linear regime, 
we may expect that the positive flux due to SDWs will aid 
and enhance the outward angular momentum transport due 
to the convective mode alone, although a further non-linear 
study is required to ascertain this. 



5 SUMMARY AND DISCUSSION 

In this paper, we have performed a detailed investigation of 
a new linear mechanism of spiral density wave excitation by 
vertical thermal convection in compressible Keplerian discs 
with superadiabatic vertical stratification, using the local 
shearing box approach. The wave excitation results from 
the velocity shear of the disc's Keplerian differential rota- 
tion. As is usually done in the shearing box, perturbations 
were decomposed into spatial Fourier harmonics, or shear- 
ing plane waves. The temporal evolution of the amplitudes of 
these waves was followed by numerical integration of the lin- 
earised hydrodynamical equations of the shearing box. Only 
non-axisymmetric perturbations were considered, as only for 
them the effects of shear are important. Three basic types of 
perturbation modes can be distinguished in the considered 
system: SDWs due mainly to compressibility, the convective 
mode due mainly to the negative vertical entropy gradient, 

4 But overall, wave perturbations and associated transport are 
not damped in a fully developed turbulence, as successive leading 
SFHs are regenerated due to non-linearity, maintaining a steady 
angular momentum flux. 



and the vortical mode due to the combined action of vertical 
stratification and the Coriolis force. However, in the present 
case of non-self-gravitating discs with superadiabatic verti- 
cal structure, the main energy-carrying mode is convection. 
In such a setup, the vortical mode can be neglected and 
hence we have concentrated primarily on the dynamics of 
SDWs and the convective mode, which both have zero po- 
tential vorticity. We first characterised the properties of the 
SDW and convective modes by deriving the dispersion re- 
lations in the presence of disc rotation and stratification, 
but neglecting shear of disc flow. In this limit of rigid ro- 
tation, the modes evolve independently, that is, once either 
of these modes has been initially excited, it does not lead 
to the excitation of another. We then demonstrated that on 
taking into account differential rotation/shear, an initially 
tightly leading SFH of the convective mode evolving in the 
disc flow, on becoming trailing, generates a corresponding 
SFH of SDWs. The main mechanism of wave generation is 
the following. Because of shear, the radial wavenumber of 
SFH varies with time, making the characteristic time-scales 
of these modes time-dependent as well. As a result, during a 
short period of time when SFHs swing from leading to trail- 
ing, the mode time-scales can become comparable to each 
other and to the shear time, thereby making an efficient en- 
ergy exchange between the SDW mode, convective mode and 
the mean disc flow possible. Consequently, SDWs are gen- 
erated in this process, mainly at the expense of shear flow 
energy, with the help of the convective mode. We also quan- 
tified the efficiency of wave generation at different azimuthal 
and vertical wavenumbers and found that it is maximal when 
these length-scales are comparable to the disc scale height. 
This, in turn, implies that SDWs excited by convection have 
weak vertical dependence, similar to SDWs generated by 
vortices (see HP09a). We calculated the angular momentum 
flux associated with non-axisymmetric SDWs generated by 
the convective mode and found that it is positive, i.e., trans- 
port is outward as opposed to that of the convective mode. 

Studied here linear coupling of SDWs and convec- 
tion is basically similar in nature to the linear coupling of 
SD Ws and vortices that has already been extensively stud- 
ied jBodo et al.ll2005l :l Mamatsashvili & Chagelishvih1l2007l . 
HP09a) and represents a special manifestation of a more 
general phenomenon of shear-induced linear coupling of per- 
turbation modes inevitably taking place in any flow with 
an inhomogeneous velocity profile, or shear. Thus, although 
convection is not generally shear-driven - in that it does 
not directly tap energy from shear but from the unstable 
superadiabatic thermal distribution - its dynamics is still 
affected by shear and shear-induced coupling to SDWs, in 
some sense, makes convection a participant in kinematic pro- 
cesses as well. 

Shearing box simulation of compressible MRI- 
turbulence in a Keplerian disc indicates that the angular 
momentum transport due to non-axisymmetric SDWs 
generated by vortical perturbations constitutes a signifi- 
cant fraction of th at associated with the total turb ulent 
Reynolds stresses ((Heinemann fc Papaloizoul l2009bl ). In 
other words, a purely hydrodynamic part of transport 
is mostly due to these SDWs and is quite appreciable. 
As mention ed in the Intro d uction , in their non-linear 
simulations, ILesur fc Qgilviel i|2010D demonstrated that, 
contrary to previous results, in fact vertical convection in 
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protoplanetary discs has a non-axisymmetric structure and 
is thus able to transport angular momentum outwards. But 
the magnitude of the corresponding a reported in their 
simulations is still small ( ,$ 10 -4 ). However, these simu- 
lations were performed in the Boussinesq/incompressible 
limit, where h igh-frequency SDWs are filtered out. In a 
related study, iKapyla, et al.l l|2010l ) solved full hydrody- 
namical equations in the shearing box with superadiabatic 
stratification without explicitly making the assumption 
of incompressibility, but the parameter regime considered 
was such as to give small Mach numbers. So, the effects of 
compressibility and, hence of SDWs, were probably small 
in their simulations and were not particularly addressed. 
As noted above, SDWs are generally known to be able to 
enhance angular momentum transport rates. For example, 
iJohnson fc Gammiel (|2005bh demonstrated that compress- 
ible simulations of vortices in 2D discs yield at least an 
order of magnitude larger transport rate, which is primarily 
due to shocks of SDWs emitted by vortices, compared 
with that in the incompressible case, which is due to 
vortices only. Based on this property, we might anticipate 
a similar situation if compressibility is taken into account 
in simulations of convectively unstable discs, where SDWs 
generated by convective motions in the non-linear regime 
could also boost outward angular momentum transport due 
to convection alone, although this needs to be further in- 
vestigated in greater detail in future numerical studies. We 
emphasise once more that the described here shear-induced 
coupling between the SDWs and convection occurs only 
when these mo des are non-a xisym metric. So, in th e earlie r 
simulations of ICabotJ l|l99cj) and IStone fc Balbusl (|l996h . 
where vertical convection had an axisymmetric structure, 
SDW generation could not be observed in spite of moderate 
Mach numbers associated with convective motions. Another 
point, which also merits investigation is where in the disc 
SDWs generated by convection pri marily cause dissipatio n 
through shock formation (see also IStone fc Balbusl 1 19961 ) . 
because in all the above-mentioned simulations convection 
was not generated self-consistently, but maintained by an 
artificially imposed heat flux. This becomes even more 
interesting in the light of the result obtained here that in 
the linear regime SDWs do not vary much with height. 
If the shocks of SDWs dissipate near the surface, this 
may stifle convection via reducing the negative entropy 
gradient, while if - near the midplane, this can then sustain 
convection. 



Finally, we would like to mention that strong vertical 
convection also occurs in FU Orionis systems during out- 
bursts as a result of hydrogen ionisation and is one of the 
important factors f or understanding the nature of the out- 
burst phenomenon (|Zhu et al.ll2009h . In this case, as shown 
by Zhu et al., convection is confined to the inner disc, ex- 
tends over the whole vertical height and causes strong den- 
sity fluctuations, because the Mach number of the convec- 
tive eddies/motions is of the order of unity. In other words, 
convection is compressible, so the question of the SDW exci- 
tation and the role of wave transport in the outburst process 
seems relevant for further investigation. 
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